はじめに

本日は地震観測データとその他の観測データを組合せて可視化する手法を学びます.

地震観測データには,アメリカ地質調査所(United States Geological Survey; USGS) Earthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)が公開するデータを用います.

USGSの地震観測データと組み合わせるその他の観測データには, テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データとスミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM(http://www.volcano.si.edu/)が公開する火山データを用います.

本日の演習では,緯度(latitude)・経度(longitude)という情報を基に複数のデータを組み合わせる手法の習得を目標とします.

プレート境界線の可視化

まず,テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データを可視化します。

プレート境界線とはプレートテクトニクス(プレート理論)により説明されるプレートとプレートの境目のことです.詳しくは以下リソースを参照してください.

http://ja.wikipedia.org/wiki/プレートテクトニクス

プレート境界線データのダウンロード

テキサス大学地球物理研究所PLATESプロジェクトのデータ公開ページ(http://www-udc.ig.utexas.edu/external/plates/data.htm)からプレート境界線データをダウンロードします.

ページの中頃に「PLATES_PlateBoundary_ArcGIS.zip」というファイルがあるので,それをダウンロードします。

Rの作業ディレクトリにディレクトリ「plate」を作成し, そのディレクトリに「PLATES_PlateBoundary_ArcGIS.zip」を展開してください.

Rで上記処理をさせる場合は以下となります.

# fileがなければダウンロードする
if(!file.exists("PLATES_PlateBoundary_ArcGIS.zip")){
  download.file(url="http://www-udc.ig.utexas.edu/external/plates/data/plate_boundaries/PLATES_PlateBoundary_ArcGIS.zip",destfile = "PLATES_PlateBoundary_ArcGIS.zip")
}
# directoryがなければファイルを展開する
if(!dir.exists("plate")){
  unzip("PLATES_PlateBoundary_ArcGIS.zip",exdir = "plate")
}

ダウンロードしたプレート境界線データはESRI Shapeフォーマットで格納されています. Shape形式はGIS(地理情報システム)で比較的よく利用されているファイルフォーマットです. 詳しくは Esri Japanのシェープファイルについて等を参照してください.

Shapeファイル処理用ライブラリのインストール

Shapeファイルを処理するため maptools ライブラリをインストールします. また,データフレーム変換用のライブラリ broomもインストールします.

#ライブラリの読み込み.なければインストールする.
if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
if(!require(maps)){
  install.packages("maps")
  library(maps)
}
if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
if(!require(broom)){
  install.packages("broom")
  library(broom)
}
#游ゴシック体を使う
if(.Platform$OS.type=="windows")
  windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
  quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))

Shapeファイルの描画

それでは早速Shapeファイルを描画しましょう.

発散型境界の描画

展開したShapeファイルを読み込み,ggplotで発散型境界を描画します.

ggplot2 は data.frame でないと処理できないので, maptoolsの処理結果(SpatialLinesDataFrame)をそのまま利用することができません. そこで,broom パッケージのtidy 関数を利用して SpatialLinesDataFrame を data.frame に変換ます.

それでは早速下記コードを実行してください.

# shapeファイルの読み込み
ridge <-readShapeSpatial("plate/ridge.shp") 
# オブジェクトのクラスを確認する
class(ridge)
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
ridge.fort<-broom::tidy(ridge)
class(ridge.fort)
## [1] "data.frame"
g <- ggplot() + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("発散型境界 (data by UTIG Plate Proj.)")
g

発散型境界線が描画されました.

世界地図と合わせて描画してみましょう.

g+borders("world")

トランスフォーム型境界の描画

つぎにトランスフォーム型境界を描画します。

transform <-readShapeSpatial("plate/transform.shp") 
transform.fort<-broom::tidy(transform)
g <- ggplot() + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("トランスフォーム型境界 (data by UTIG Plate Proj.)")
g

これも世界地図と合わせて描画してみましょう.

g+borders("world")

収束型境界の描画

つぎに収束型境界を描画します。

trench <-readShapeSpatial("plate/trench.shp") 
trench.fort<-broom::tidy(trench)
g <- ggplot() + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("収束型境界 (data by UTIG Plate Proj.)")
g

収束型境界についても,世界地図と合わせて描画してみましょう.

g+borders("world")

すべての境界線の描画

それでは発散型,トランスフォーム型,収束型境界線を合わせて描画します.

ridge <-readShapeSpatial("plate/ridge.shp") 
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp") 
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp") 
trench.fort<-broom::tidy(trench)
g <- ggplot()
g <- g + theme_bw(base_family = "yugo")+ggtitle("プレート境界 (data by UTIG Plate Proj.)")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
g

世界地図と合わせます.

g + borders(database = "world")

地震はプレート境界線で生じる歪が原因となり発生するといわれています.

日本で地震が多い理由も納得できますね.

地震観測データとプレート境界線の組合せ

地震観測データの描画

アメリカ地質調査所(United States Geological Survey; USGS) のEarthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)の地震観測データとテキサス大学地球物理研究所PLATESプロジェクトのプレート境界データを合成します.

まずは地震観測データを描画します. データは前回の可視化で用いた「20110311JSTquery.csv」を使います.

csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
#地震をプロット(色で深度,サイズでマグニチュードを表現)
g<-ggplot(csv,aes(x=longitude,y=latitude,color=depth,size=mag))
g<-g+geom_point()
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ(data by USGS Earthquake Hazards Prog. Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く,深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
g

ここまでは前回と同じですね.

地震観測データとプレート境界線データの描画

それでは地震観測データとプレート境界線データと合わせて可視化してみましょう.

ポイントはこれまで描画データを指定していたggplot()関数では何も指定せずに, 点や線を描画する geom_point や geom_path で描画データを指定している点です. こうすることで複数のデータを用いた可視化が可能となります(ggplot 関数でデータを指定すると全体でデータが共有されることになる).

# 地震観測データの読み込み
csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
# プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp") 
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp") 
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp") 
trench.fort<-broom::tidy(trench)

# 描画データを指定しない
g<-ggplot()
# geom_point で描画するデータ(地震観測データ)を指定している.
g<-g+geom_point(data = csv,mapping = aes(x=longitude,y=latitude,color=depth,size=mag))
# 世界地図の描画
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ(data by USGS Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く,深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
# タイトルを上書き
g <- g +ggtitle("プレート境界と地震データ(data by UTIG Plate Proj. and USGS)")
#発散型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
g

プレート境界線近辺で地震が発生していることがわかります.

地震発生メカニズム(気象庁,地震発生のしくみ)を実際のデータで確認することができました.

プレート境界線データと火山データの組合せ

つぎにプレート境界線データと火山データを組合せて可視化してみます.

火山もプレートテクトニクスに関連が強い地形です.(参考:Wikipdia,火山

火山データには,スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM(http://www.volcano.si.edu/)が公開する火山データを用います.

火山データのダウンロードとフォーマット変換

火山データのダウンロードするため,スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMサイト(http://www.volcano.si.edu/)へアクセスします.

「Database」 > 「HOLOCENE SPREADSHEET」をクリックします.

「GVP_Volcano_List.xls」というファイル名でダウンロードが始まります.

ダウンロードが完了したら,「GVP_Volcano_List.xls」をCSVフォーマットに変換します.

フォーマット変換にはExcelを利用します.Excelでファイルを開き,「名前を付けて保存」でCSVフォーマットとして保存します.

変換後のCSVファイル「GVP_Volcano_List.csv」をRの作業ディレクトリに移動します.

火山データの描画

それではさっそく火山データを読み込んでみましょう.

#データを読み込む
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1,stringsAsFactors = FALSE)
#構造を確認する
str(volc)
## 'data.frame':    1520 obs. of  13 variables:
##  $ Volcano.Number      : int  210010 210020 210030 210040 211001 211003 211004 211010 211020 211030 ...
##  $ Volcano.Name        : chr  "West Eifel Volcanic Field" "Chaine des Puys" "Olot Volcanic Field" "Calatrava Volcanic Field" ...
##  $ Country             : chr  "Germany" "France" "Spain" "Spain" ...
##  $ Primary.Volcano.Type: chr  "Maar(s)" "Lava dome(s)" "Pyroclastic cone(s)" "Pyroclastic cone(s)" ...
##  $ Activity.Evidence   : chr  "Eruption Dated" "Eruption Dated" "Evidence Credible" "Eruption Dated" ...
##  $ Last.Known.Eruption : chr  "8300 BCE" "4040 BCE" "Unknown" "3600 BCE" ...
##  $ Region              : chr  "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" ...
##  $ Subregion           : chr  "Western Europe" "Western Europe" "Western Europe" "Western Europe" ...
##  $ Latitude            : num  50.2 45.8 42.2 38.9 43.2 ...
##  $ Longitude           : num  6.85 2.97 2.53 -4.02 10.87 ...
##  $ Elevation..m.       : int  600 1464 893 1117 500 800 949 458 1281 789 ...
##  $ Dominant.Rock.Type  : chr  "Foidite" "Basalt / Picro-Basalt" "Trachybasalt / Tephrite Basanite" "Basalt / Picro-Basalt" ...
##  $ Tectonic.Setting    : chr  "Rift zone / Continental crust (>25 km)" "Rift zone / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" ...
#要約を確認する
summary(volc)
##  Volcano.Number   Volcano.Name         Country         
##  Min.   :210010   Length:1520        Length:1520       
##  1st Qu.:261148   Class :character   Class :character  
##  Median :300020   Mode  :character   Mode  :character  
##  Mean   :296768                                        
##  3rd Qu.:342123                                        
##  Max.   :600000                                        
##  Primary.Volcano.Type Activity.Evidence  Last.Known.Eruption
##  Length:1520          Length:1520        Length:1520        
##  Class :character     Class :character   Class :character   
##  Mode  :character     Mode  :character   Mode  :character   
##                                                             
##                                                             
##                                                             
##     Region           Subregion            Latitude        Longitude      
##  Length:1520        Length:1520        Min.   :-78.50   Min.   :-179.97  
##  Class :character   Class :character   1st Qu.: -7.13   1st Qu.: -78.15  
##  Mode  :character   Mode  :character   Median : 13.84   Median :  38.39  
##                                        Mean   : 14.05   Mean   :  23.85  
##                                        3rd Qu.: 42.00   3rd Qu.: 139.04  
##                                        Max.   : 85.61   Max.   : 179.58  
##  Elevation..m.     Dominant.Rock.Type Tectonic.Setting  
##  Min.   :-4200.0   Length:1520        Length:1520       
##  1st Qu.:  680.8   Class :character   Class :character  
##  Median : 1462.0   Mode  :character   Mode  :character  
##  Mean   : 1677.0                                        
##  3rd Qu.: 2344.8                                        
##  Max.   : 6879.0
#データを確認する
head(volc)
##   Volcano.Number              Volcano.Name Country Primary.Volcano.Type
## 1         210010 West Eifel Volcanic Field Germany              Maar(s)
## 2         210020           Chaine des Puys  France         Lava dome(s)
## 3         210030       Olot Volcanic Field   Spain  Pyroclastic cone(s)
## 4         210040  Calatrava Volcanic Field   Spain  Pyroclastic cone(s)
## 5         211001                Larderello   Italy  Explosion crater(s)
## 6         211003                   Vulsini   Italy              Caldera
##   Activity.Evidence Last.Known.Eruption                         Region
## 1    Eruption Dated            8300 BCE Mediterranean and Western Asia
## 2    Eruption Dated            4040 BCE Mediterranean and Western Asia
## 3 Evidence Credible             Unknown Mediterranean and Western Asia
## 4    Eruption Dated            3600 BCE Mediterranean and Western Asia
## 5 Eruption Observed             1282 CE Mediterranean and Western Asia
## 6 Eruption Observed             104 BCE Mediterranean and Western Asia
##        Subregion Latitude Longitude Elevation..m.
## 1 Western Europe   50.170      6.85           600
## 2 Western Europe   45.775      2.97          1464
## 3 Western Europe   42.170      2.53           893
## 4 Western Europe   38.870     -4.02          1117
## 5          Italy   43.250     10.87           500
## 6          Italy   42.600     11.93           800
##                 Dominant.Rock.Type
## 1                          Foidite
## 2            Basalt / Picro-Basalt
## 3 Trachybasalt / Tephrite Basanite
## 4            Basalt / Picro-Basalt
## 5                No Data (checked)
## 6          Trachyte / Trachydacite
##                               Tectonic.Setting
## 1       Rift zone / Continental crust (>25 km)
## 2       Rift zone / Continental crust (>25 km)
## 3      Intraplate / Continental crust (>25 km)
## 4      Intraplate / Continental crust (>25 km)
## 5 Subduction zone / Continental crust (>25 km)
## 6 Subduction zone / Continental crust (>25 km)
#富士山のデータを確認する
volc[volc$Volcano.Name=="Fujisan",]
##     Volcano.Number Volcano.Name Country Primary.Volcano.Type
## 598         283030      Fujisan   Japan        Stratovolcano
##     Activity.Evidence Last.Known.Eruption                  Region
## 598 Eruption Observed             1708 CE Japan, Taiwan, Marianas
##     Subregion Latitude Longitude Elevation..m.    Dominant.Rock.Type
## 598    Honshu   35.361   138.728          3776 Basalt / Picro-Basalt
##                                 Tectonic.Setting
## 598 Subduction zone / Continental crust (>25 km)

火山データの位置が緯度(Latitude),経度(Longitude)に記録されていることがわかります.

それでは火山データを可視化してみましょう.ここでは散布図 geom_point を用いていますが,単なる点だと火山ぽくないので,「shape=24」を指定することで点を三角に変更しています.

#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")

geom_point で指定可能な shape には他にも以下のような形があります.

g<-ggplot(data=data.frame(x=0:24))+scale_shape_identity()
g<-g+geom_point(mapping=aes(x=x,y=0,shape=x),size=5)
g<-g+geom_point(mapping=aes(x=x,y=1,shape=x),size=5,fill="red")
g<-g+xlab("shape types")
g

今回はshapeの24番(色塗り可能な三角形)を指定することで火山を表現しています.

スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMの火山データにはその他にも様々な情報が記録されています.例えば,「Elevation..m.」には火山の標高が記録されています.

可視化してみると次のようになります.

volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude,color=Elevation..m.,fill=Elevation..m.),shape=24)+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program)")

火山データとプレート境界線データの描画

それでは火山データとプレート境界線データを合わせて描画してみましょう.

まず,火山データと世界地図を描画します.

volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
g<-ggplot()
g<-g+borders(database = "world")
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
g

この描画結果にプレート境界線を追加すればよいわけです.

地震観測データとプレート境界線データの描画を参考にすると, 下記のようにコードを書くことができると思います.要は描画単位(geom_pointやgeom_pathごと)に可視化するデータを指定すればよいのです.

#プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp") 
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp") 
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp") 
trench.fort<-broom::tidy(trench)

#火山データの読み込み
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)

g<-ggplot()
#地図の描画
g<-g+borders(database = "world")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
#火山の描画
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
g

このように火山もプレート境界線近辺に多くできていることがわかります.

まとめ

本日は地球観測データと科学データを複数組合せて可視化する手法を学びました. 具体的には,緯度(latitude)・経度(longitude)をキーに,地震観測データ,プレート境界線データ,火山データを関連付けて可視化しました.このように複数のデータを組合せて可視化するにはデータを関連付けるためのキーが必要となります.今回の演習で感じてもらえたと思いますが,緯度・経度は複数データを組み合わせる際のキー候補としやすい情報の一つです.

次回は今回の内容に加えて地理空間情報の可視化に取り組みます.

課題

 以下の問題を解いてください. その際に用いたRスクリプトファイル,結果として得られた画像ファイルをZIPでまとめて提出してください.

※Rスクリプトファイルは可視化毎に作成してください.可視化に必要な内容のみを含めるようにしてください. ※データファイルは不要です.

【問題】地震,プレート境界線,火山データのすべてを可視化してください

ヒント:

  • 火山データの標高は可視化しなくて構いません
  • 見やすくなるよう描画順や色等を工夫してください

可視化例:

【発展問題】上記内容を特定の国に絞って可視化してください

ヒント:

  • ggplot2 の可視化データはすべてdata.frameです
  • subsetを使います

例:

戻る


クリエイティブ・コモンズ・ライセンス
Masaharu Hayashi を著作者とするこの 作品 は クリエイティブ・コモンズの 表示 4.0 国際 ライセンスで提供されています。